<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Varying Parameters analysis</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S2 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S3 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S4 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: pre-wrap;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: pre;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: pre;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S5 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S6 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S7 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S8 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S12 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S13 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: pre-wrap;  }
.variableValue { width: 100% !important; }
.embeddedOutputsMatrixElement,.eoOutputWrapper .matrixElement {    min-height: 18px;    box-sizing: border-box;}
.embeddedOutputsMatrixElement .matrixElement,.eoOutputWrapper  .matrixElement,.rtcDataTipElement .matrixElement {    position: relative;}
.matrixElement .variableValue,.rtcDataTipElement .matrixElement .variableValue {    white-space: pre;    display: inline-block;    vertical-align: top;    overflow: hidden;}
.embeddedOutputsMatrixElement.inlineElement {}
.embeddedOutputsMatrixElement.inlineElement .topHeaderWrapper {    display: none;}
.embeddedOutputsMatrixElement.inlineElement .veTable .body {    padding-top: 0 !important;    max-height: 100px;}
.inlineElement .matrixElement {    max-height: 300px;}
.embeddedOutputsMatrixElement.rightPaneElement {}
.rightPaneElement .matrixElement,.rtcDataTipElement .matrixElement {    overflow: hidden;    padding-left: 9px;}
.rightPaneElement .matrixElement {    margin-bottom: -1px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer,.eoOutputWrapper .matrixElement .valueContainer,.rtcDataTipElement .matrixElement .valueContainer {    white-space: nowrap;    margin-bottom: 3px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis.hide,.embeddedOutputsMatrixElement .matrixElement .verticalEllipsis.hide,.eoOutputWrapper .matrixElement .valueContainer .horizontalEllipsis.hide,.eoOutputWrapper .matrixElement .verticalEllipsis.hide,.rtcDataTipElement .matrixElement .valueContainer .horizontalEllipsis.hide,.rtcDataTipElement .matrixElement .verticalEllipsis.hide {    display: none;}
.embeddedOutputsVariableMatrixElement .matrixElement .valueContainer.hideEllipses .verticalEllipsis, .embeddedOutputsVariableMatrixElement .matrixElement .valueContainer.hideEllipses .horizontalEllipsis {    display:none;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis,.eoOutputWrapper .matrixElement .valueContainer .horizontalEllipsis {    margin-bottom: -3px;}
.eoOutputWrapper .embeddedOutputsVariableMatrixElement .matrixElement .valueContainer {    cursor: default !important;}
/* * Ellipses as base64 for HTML export. */.matrixElement .horizontalEllipsis,.rtcDataTipElement .matrixElement .horizontalEllipsis {    display: inline-block;    margin-top: 3px;    /* base64 encoded version of images-liveeditor/HEllipsis.png */    width: 30px;    height: 12px;    background-repeat: no-repeat;    background-image: url("");}
.matrixElement .verticalEllipsis,.textElement .verticalEllipsis,.rtcDataTipElement .matrixElement .verticalEllipsis,.rtcDataTipElement .textElement .verticalEllipsis {    margin-left: 35px;    /* base64 encoded version of images-liveeditor/VEllipsis.png */    width: 12px;    height: 30px;    background-repeat: no-repeat;    background-image: url("");}</style></head><body><div class = rtcContent><h1  class = 'S0'><span>Varying Parameters analysis</span></h1><div  class = 'S1'><span style=' font-weight: bold;'>Author(s): Vanja Vlasov, Systems Biochemistry Group, LCSB, University of Luxembourg.</span></div><div  class = 'S1'><span style=' font-weight: bold;'>Reviewer(s): Thomas Pfau, Systems Biology Group, LSRU, University of Luxembourg.</span></div><div  class = 'S1'><span style=' font-weight: bold;'>Ines Thiele, Molecular Systems Physiology Group, LCSB, University of Luxembourg</span></div><div  class = 'S1'><span style=' font-weight: bold;'>Almut Heinken, Molecular Systems Physiology Group, LCSB, University of Luxembourg</span></div><div  class = 'S1'><span style=' font-weight: bold;'></span></div><div  class = 'S1'><span>In this tutorial, we show how computations are performed by varying one or two parameters over a fixed range of numerical values.</span></div><h2  class = 'S2'><span>EQUIPMENT SETUP</span></h2><div  class = 'S1'><span>If necessary, initialise the cobra toolbox:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S3'><span style="white-space: pre"><span >initCobraToolbox(false) </span><span style="color: rgb(2, 128, 9);">% false, as we don't want to update</span></span></div><div  class = 'S4'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="57695B64" data-testid="output_0" data-width="420" data-height="899" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">      _____   _____   _____   _____     _____     |
     /  ___| /  _  \ |  _  \ |  _  \   / ___ \    |   COnstraint-Based Reconstruction and Analysis
     | |     | | | | | |_| | | |_| |  | |___| |   |   The COBRA Toolbox - 2017
     | |     | | | | |  _  { |  _  /  |  ___  |   |
     | |___  | |_| | | |_| | | | \ \  | |   | |   |   Documentation:
     \_____| \_____/ |_____/ |_|  \_\ |_|   |_|   |   <a href="http://opencobra.github.io/cobratoolbox" style="white-space: pre; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">http://opencobra.github.io/cobratoolbox</a>
                                                  | 

 &gt; Checking if git is installed ...  Done.
 &gt; Checking if the repository is tracked using git ...  Done.
 &gt; Checking if curl is installed ...  Done.
 &gt; Checking if remote can be reached ...  Done.
 &gt; Initializing and updating submodules ... Done.
 &gt; Adding all the files of The COBRA Toolbox ...  Done.
 &gt; Define CB map output... set to svg.
 &gt; Retrieving models ...   Done.
 &gt; TranslateSBML is installed and working properly.
 &gt; Configuring solver environment variables ...
   - [*---] ILOG_CPLEX_PATH: C:\Program Files\IBM\ILOG\CPLEX_Studio1263\cplex\matlab\x64_win64
   - [*---] GUROBI_PATH: C:\gurobi650\win64\matlab
   - [*---] TOMLAB_PATH: C:\tomlab\
   - [----] MOSEK_PATH :  --&gt; set this path manually after installing the solver ( see <a href="https://opencobra.github.io/cobratoolbox/docs/solvers.html" style="white-space: pre; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">instructions</a> )
   Done.
 &gt; Checking available solvers and solver interfaces ... Done.
 &gt; Setting default solvers ... Done.
 &gt; Saving the MATLAB path ... Done.
   - The MATLAB path was saved in the default location.

 &gt; Summary of available solvers and solver interfaces

					Support           LP 	 MILP 	   QP 	 MIQP 	  NLP
	----------------------------------------------------------------------
	cplex_direct 	full          	    0 	    0 	    0 	    0 	    -
	dqqMinos     	full          	    0 	    - 	    - 	    - 	    -
	glpk         	full          	    1 	    1 	    - 	    - 	    -
	gurobi       	full          	    1 	    1 	    1 	    1 	    -
	ibm_cplex    	full          	    0 	    0 	    0 	    - 	    -
	matlab       	full          	    1 	    - 	    - 	    - 	    1
	mosek        	full          	    0 	    0 	    0 	    - 	    -
	pdco         	full          	    1 	    - 	    1 	    - 	    -
	quadMinos    	full          	    0 	    - 	    - 	    - 	    0
	tomlab_cplex 	full          	    1 	    1 	    1 	    1 	    -
	qpng         	experimental  	    - 	    - 	    1 	    - 	    -
	tomlab_snopt 	experimental  	    - 	    - 	    - 	    - 	    1
	gurobi_mex   	legacy        	    0 	    0 	    0 	    0 	    -
	lindo_old    	legacy        	    0 	    - 	    - 	    - 	    -
	lindo_legacy 	legacy        	    0 	    - 	    - 	    - 	    -
	lp_solve     	legacy        	    1 	    - 	    - 	    - 	    -
	opti         	legacy        	    0 	    0 	    0 	    0 	    0
	----------------------------------------------------------------------
	Total        	-             	    6 	    3 	    4 	    2 	    2

 + Legend: - = not applicable, 0 = solver not compatible or not installed, 1 = solver installed.


 &gt; You can solve LP problems using: 'glpk' - 'gurobi' - 'matlab' - 'pdco' - 'tomlab_cplex' - 'lp_solve' 
 &gt; You can solve MILP problems using: 'glpk' - 'gurobi' - 'tomlab_cplex' 
 &gt; You can solve QP problems using: 'gurobi' - 'pdco' - 'tomlab_cplex' - 'qpng' 
 &gt; You can solve MIQP problems using: 'gurobi' - 'tomlab_cplex' 
 &gt; You can solve NLP problems using: 'matlab' - 'tomlab_snopt' 

 &gt; Checking for available updates ...
 --&gt; You cannot update your fork using updateCobraToolbox(). [535a88 @ develop].
     Please use the MATLAB.devTools (<a href="https://github.com/opencobra/MATLAB.devTools" style="white-space: pre; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">https://github.com/opencobra/MATLAB.devTools</a>).</div></div></div></div></div><div  class = 'S5'><span>For solving line</span><span>ar programming problems in the analysis, certain solvers are required:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >changeCobraSolver (</span><span style="color: rgb(170, 4, 249);">'gurobi'</span><span >, </span><span style="color: rgb(170, 4, 249);">'all'</span><span >, 1);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%changeCobraSolver ('glpk', 'all', 1);</span></span></div><div  class = 'S4'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="9837BD7C" data-testid="output_1" data-width="420" data-height="73" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> &gt; Solver for LPproblems has been set to glpk.
 &gt; Solver for MILPproblems has been set to glpk.
 &gt; Solver glpk not supported for problems of type MIQP. Currently used: tomlab_cplex 
 &gt; Solver glpk not supported for problems of type NLP. Currently used: matlab 
 &gt; Solver glpk not supported for problems of type QP. Currently used: qpng </div></div></div></div></div><div  class = 'S5'><span>The present tutoria</span><span>l can ru</span><span>n </span><span>with </span><span style=' font-family: monospace;'>'glpk</span><span style=' font-family: monospace;'>'</span><span> package, which does not require additional installation and configuration. Although, for the analysis of large models is recommended to use the </span><span style=' font-family: monospace;'>'gurobi'</span><span> package.</span></div><h2  class = 'S2'><span>PROCEDURE</span></h2><div  class = 'S1'><span>Before proceeding with the simulations, the path for the model needs to be set up. In thi</span><span>s tutorial, the used</span><span> model is the generic model of human metabolism, Recon 3 [1]. Therefore, we assume, that the cellular objectives include energy production or optimisation of uptake rates and by-product secretion for various physiological functions of the human body. If Recon 3 is not available, please use Recon 2.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%For Recon3D Change the model</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelFileName = </span><span style="color: rgb(170, 4, 249);">'Recon2.0model.mat'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelDirectory = getDistributedModelFolder(modelFileName); </span><span style="color: rgb(2, 128, 9);">%Look up the folder for the distributed Models.</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelFileName= [modelDirectory filesep modelFileName]; </span><span style="color: rgb(2, 128, 9);">% Get the full path. Necessary to be sure, that the right model is loaded</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >model = readCbModel(modelFileName);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'></div></div></div><div  class = 'S5'><span>If Recon2 is used, the reaction nomenclature needs to be adjusted.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >model.rxns(find(ismember(model.rxns,</span><span style="color: rgb(170, 4, 249);">'EX_glc(e)'</span><span >)))={</span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >model.rxns(find(ismember(model.rxns,</span><span style="color: rgb(170, 4, 249);">'EX_o2(e)'</span><span >)))={</span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >};</span></span></div></div></div><h2  class = 'S2'><span>TROUBLESHOOTING</span></h2><div  class = 'S1'><span>If there are multiple energy sources available in the model; Specifying more constraints is necessary. If we do not do that, we will have additional carbon and oxygen energy sources available in the cell and the maximal ATP production. </span></div><div  class = 'S1'><span>To avoid this issue, all external carbon sources need to be closed.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%Closing the uptake of all energy and oxygen sources</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >i=1:length(model.rxns)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >strncmp(model.rxns{i},</span><span style="color: rgb(170, 4, 249);">'EX_'</span><span >,3)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >        model.subSystems{i}=</span><span style="color: rgb(170, 4, 249);">'Exchange/demand reaction'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >idx=strmatch(</span><span style="color: rgb(170, 4, 249);">'Exchange/demand reaction'</span><span >, model.subSystems);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >c=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >i=1:length(idx)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >model.lb(idx(i))~=0</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >        c=c+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >        uptakes{c}=model.rxns{idx(i)};</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelalter = model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelalter = changeRxnBounds(modelalter, uptakes, 0, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% The alternative way to do that, in case you were using another large model, </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% that does not contain defined Subsystem is</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% to find uptake exchange reactions with following codes:</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% [selExc, selUpt] = findExcRxns(model);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% uptakes = model.rxns(selUpt);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Selecting from the exchange uptake reactions those </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% which contain at least 1 carbon in the metabolites included in the reaction:</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% subuptakeModel = extractSubNetwork(model, uptakes);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% hiCarbonRxns = findCarbonRxns(subuptakeModel,1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Closing the uptake of all the carbon sources</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelalter = model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelalter = changeRxnBounds(modelalter, hiCarbonRxns, 0, 'l');</span></span></div></div></div><h2  class = 'S2'><span>Robustness analysis</span></h2><div  class = 'S1'><span>Robustness analysis is applied to estimate and visualise how changes in the concentration of an environmental parameter (exchange rate) or internal reaction effect on the objective [2]. If we are interested in varying </span><span texencoding="$v_{j}" style="vertical-align:-6px"><img src="" width="14.5" height="20" /></span><span> between two values, i.e., </span><span texencoding="$v_{j,min}" style="vertical-align:-6px"><img src="" width="32" height="20" /></span><span> and </span><span texencoding="$v_{j,max}" style="vertical-align:-6px"><img src="" width="34" height="20" /></span><span>, we can solve </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">l</span><span> optimisation problems:</span></div><div  class = 'S10'><span texencoding="\begin{array}{ll}
\max Z_{k}= c^{T}v \\
\text{s.t.} &amp; k =1,...,l,\\
&amp; Sv=0,\\
\text{fixing} &amp; v_{j}=v_{j,min}+ \frac{(k-1)}{(l-1)}*(v_{j,max}-v_{j,min})\\
\text{constraints} &amp; v_{i,min}\leq v_{i} \leq v_{i,max}, i=1,...,n,i \neq j
 \end{array}" style="vertical-align:-58px"><img src="" width="317.5" height="127" /></span></div><div  class = 'S1'><span></span></div><div  class = 'S1'><span>The function </span><span style=' font-family: monospace;'>robustnessAnalysis()</span><span> is used for this analysis:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% [controlFlux, objFlux] = robustnessAnalysis(model, controlRxn, nPoints,...</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%      plotResFlag, objRxn,objType)</span></span></div></div></div><div  class = 'S5'><span>where inputs are a COBRA model, a reaction that has been analysed and optional inputs:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% INPUTS</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% model         COBRA model structure</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% controlRxn    Reaction of interest whose value is to be controlled</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% OPTIONAL INPUTS</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% nPoints       Number of points to show on plot (Default = 20)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% plotResFlag   Plot results (Default true)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% objRxn        Objective reaction to be maximized</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%               (Default = whatever is defined in model)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% objType       Maximize ('max') or minimize ('min') objective</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%               (Default = 'max')</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% OUTPUTS</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% controlFlux   Flux values within the range of the maximum and minimum for</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%               a reaction of interest</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% objFlux       Optimal values of objective reaction at each control</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%               reaction flux value</span></span></div></div></div><div  class = 'S5'><span>Here, we will investigate how robust the maximal ATP production of the network (</span><span>i.e., the maximal flux through '</span><span style=' font-family: monospace;'>DM_atp_c_'</span><span>) is with respect to varying glucose uptake rates and fixed oxygen uptake. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >modelrobust = modelalter;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelrobust = changeRxnBounds(modelrobust, </span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >, -17, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >AtpRates = zeros(21, 1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >i = 0:20</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    modelrobust = changeRxnBounds(modelrobust, </span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >, -i, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    modelrobust = changeObjective(modelrobust, </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    FBArobust = optimizeCbModel(modelrobust, </span><span style="color: rgb(170, 4, 249);">'max'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    AtpRates(i+1) = FBArobust.f;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre"><span >plot (1:21, AtpRates)</span></span></div><div  class = 'S4'><div class="inlineElement eoOutputWrapper embeddedOutputsWarningElement" uid="49E65269" data-testid="output_2" data-width="420" data-height="44" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType" style="white-space: normal; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"><div class="diagnosticMessage-messagePart" style="white-space: pre-wrap; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;">Warning: MATLAB has disabled some advanced graphics rendering features by switching to software OpenGL. For more information, click here.</div><div class="diagnosticMessage-stackPart" style="white-space: pre; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"></div></div></div></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >xlabel(</span><span style="color: rgb(170, 4, 249);">'Flux through EX-glc-D[e]'</span><span >)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre"><span >ylabel(</span><span style="color: rgb(170, 4, 249);">'Objective function'</span><span >)</span></span></div><div  class = 'S4'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="21E9088F" data-testid="output_3" style="width: 450px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><div  class = 'S5'><span>We can also investigate the robustness of the maximal ATP production when the available glucose amount is fixed, while different levels of oxygen are available.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >modelrobustoxy = modelalter;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelrobustoxy = changeRxnBounds(modelrobustoxy, </span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >, -20, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >AtpRatesoxy = zeros(21, 1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >i = 0:20</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    modelrobustoxy = changeRxnBounds(modelrobustoxy, </span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >, -i, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    modelrobustoxy = changeObjective(modelrobustoxy, </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    FBArobustoxy = optimizeCbModel(modelrobustoxy, </span><span style="color: rgb(170, 4, 249);">'max'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    AtpRatesoxy(i+1) = FBArobustoxy.f;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >plot (1:21, AtpRatesoxy)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >xlabel(</span><span style="color: rgb(170, 4, 249);">'Flux through EX-o2[e]'</span><span >)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre"><span >ylabel(</span><span style="color: rgb(170, 4, 249);">'Objective function'</span><span >)</span></span></div><div  class = 'S4'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="468B79E6" data-testid="output_4" style="width: 450px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><ul  class = 'S12'><li  class = 'S13'><span>   </span><span style=' font-weight: bold;'>    </span><span style=' font-weight: bold;'>Double robust analysis</span></li></ul><div  class = 'S1'><span>Performs robustness analysis for a pair of reactions of interest and an objective of interest. The double robust analysis is implemented with the function </span><span style=' font-family: monospace;'>doubleRobustnessAnalysis()</span><span>.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% [controlFlux1, controlFlux2, objFlux] = doubleRobustnessAnalysis(model,...</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%  controlRxn1, controlRxn2, nPoints, plotResFlag, objRxn, objType)</span></span></div></div></div><div  class = 'S5'><span>The inputs are a COBRA model, two reactions for the analysis and optional inputs:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%INPUTS</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% model         COBRA model to analyse,</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% controlRxn1   The first reaction for the analysis,</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% controlRxn2   The second reaction for the analysis;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%OPTIONAL INPUTS</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% nPoints       The number of flux values per dimension (Default = 20)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% plotResFlag   Indicates whether the result should be plotted (Default = true)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% objRxn        is objective to be used in the analysis (Default = whatever</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%               is defined in model)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% objType       Direction of the objective (min or max)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%               (Default = 'max')</span></span></div></div></div><div  class = 'S5'><span></span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >modeldrobustoxy = modelalter;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modeldrobustoxy = changeRxnBounds(modeldrobustoxy, </span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >, -20, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modeldrobustoxy = changeRxnBounds(modeldrobustoxy, </span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >, -17, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >[controlFlux1, controlFlux2, objFlux] = doubleRobustnessAnalysis(modeldrobustoxy,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >, </span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >, 10, 1, </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >, </span><span style="color: rgb(170, 4, 249);">'max'</span><span >)</span></span></div><div  class = 'S4'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="CE80CF32" data-testid="output_5" data-width="420" data-height="31" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Double robustness analysis in progress ...
1%      [                                        ]2%      [                                        ]3%      [.                                       ]4%      [.                                       ]5%      [..                                      ]6%      [..                                      ]7%      [..                                      ]8%      [...                                     ]9%      [...                                     ]10%     [....                                    ]11%     [....                                    ]12%     [....                                    ]13%     [.....                                   ]14%     [.....                                   ]15%     [......                                  ]16%     [......                                  ]17%     [......                                  ]18%     [.......                                 ]19%     [.......                                 ]20%     [........                                ]21%     [........                                ]22%     [........                                ]23%     [.........                               ]24%     [.........                               ]25%     [..........                              ]26%     [..........                              ]27%     [..........                              ]28%     [...........                             ]28%     [...........                             ]30%     [............                            ]31%     [............                            ]32%     [............                            ]33%     [.............                           ]34%     [.............                           ]35%     [..............                          ]36%     [..............                          ]37%     [..............                          ]38%     [...............                         ]39%     [...............                         ]40%     [................                        ]41%     [................                        ]42%     [................                        ]43%     [.................                       ]44%     [.................                       ]45%     [..................                      ]46%     [..................                      ]47%     [..................                      ]48%     [...................                     ]49%     [...................                     ]50%     [....................                    ]51%     [....................                    ]52%     [....................                    ]53%     [.....................                   ]54%     [.....................                   ]55%     [......................                  ]56%     [......................                  ]56%     [......................                  ]57%     [......................                  ]59%     [.......................                 ]60%     [........................                ]61%     [........................                ]62%     [........................                ]63%     [.........................               ]64%     [.........................               ]65%     [..........................              ]66%     [..........................              ]67%     [..........................              ]68%     [...........................             ]69%     [...........................             ]70%     [............................            ]71%     [............................            ]72%     [............................            ]73%     [.............................           ]74%     [.............................           ]75%     [..............................          ]76%     [..............................          ]77%     [..............................          ]78%     [...............................         ]79%     [...............................         ]80%     [................................        ]81%     [................................        ]82%     [................................        ]83%     [.................................       ]84%     [.................................       ]85%     [..................................      ]86%     [..................................      ]87%     [..................................      ]88%     [...................................     ]89%     [...................................     ]90%     [....................................    ]91%     [....................................    ]92%     [....................................    ]93%     [.....................................   ]94%     [.....................................   ]95%     [......................................  ]96%     [......................................  ]97%     [......................................  ]98%     [....................................... ]99%     [....................................... ]100%    [........................................]</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsMatrixElement" uid="30B78F01" data-testid="output_6" data-width="420" style="width: 450px;"><div class="matrixElement veSpecifier"><div class="veVariableName variableNameElement" style="width: 420px;"><span>controlFlux1 = </span><span class="veVariableValueSummary"></span></div><div class="valueContainer" data-layout="{&quot;columnWidth&quot;:72,&quot;totalColumns&quot;:&quot;1&quot;,&quot;totalRows&quot;:&quot;10&quot;,&quot;charsPerColumn&quot;:10}"><div class="variableValue" style="width: 74px;">  -20.0000
  -17.7225
  -15.4451
  -13.1676
  -10.8902
   -8.6127
   -6.3353
   -4.0578
   -1.7804
    0.4971
</div><div class="horizontalEllipsis hide"></div><div class="verticalEllipsis hide"></div></div></div></div><div class="inlineElement eoOutputWrapper embeddedOutputsMatrixElement" uid="E00A069A" data-testid="output_7" data-width="420" style="width: 450px;"><div class="matrixElement veSpecifier"><div class="veVariableName variableNameElement" style="width: 420px;"><span>controlFlux2 = </span><span class="veVariableValueSummary"></span></div><div class="valueContainer" data-layout="{&quot;columnWidth&quot;:72,&quot;totalColumns&quot;:&quot;1&quot;,&quot;totalRows&quot;:&quot;10&quot;,&quot;charsPerColumn&quot;:10}"><div class="variableValue" style="width: 74px;">  -17.0000
  -15.1111
  -13.2222
  -11.3333
   -9.4444
   -7.5556
   -5.6667
   -3.7778
   -1.8889
    0.0000
</div><div class="horizontalEllipsis hide"></div><div class="verticalEllipsis hide"></div></div></div></div><div class="inlineElement eoOutputWrapper embeddedOutputsMatrixElement" uid="E3211CAA" data-testid="output_8" data-width="420" style="width: 450px;"><div class="matrixElement veSpecifier"><div class="veVariableName variableNameElement" style="width: 420px;"><span>objFlux = </span><span class="veVariableValueSummary"></span></div><div class="valueContainer" data-layout="{&quot;columnWidth&quot;:72,&quot;totalColumns&quot;:&quot;10&quot;,&quot;totalRows&quot;:&quot;10&quot;,&quot;charsPerColumn&quot;:10}"><div class="variableValue" style="width: 362px;">  126.2944  116.7061  107.1179   97.5296   87.9413   78.3531   68.7648   59.1765   49.5883   40.0000
  121.7395  112.1512  102.5630   92.9747   83.3864   73.7982   64.2099   54.6216   45.0334   35.4451
  117.1846  107.5963   98.0081   88.4198   78.8315   69.2433   59.6550   50.0667   40.4785   30.8902
  112.6297  103.0414   93.4532   83.8649   74.2766   64.6884   55.1001   45.5118   35.9236   26.3353
  108.0748   98.4865   88.8983   79.3100   69.7217   60.1335   50.5452   40.9569   31.3686   21.7804
  103.5199   93.9316   84.3434   74.7551   65.1668   55.5785   45.9903   36.4020   26.8137   17.2255
   98.9650   89.3767   79.7884   70.2002   60.6119   51.0236   41.4354   31.8471   22.2588   12.6706
   94.4101   84.8218   75.2335   65.6453   56.0570   46.4687   36.8805   27.2922   17.7039    8.1157
   87.7466   78.7732   69.7997   60.8263   51.5021   41.9138   32.3256   22.7373   13.1490    3.5608
   33.5029         0         0         0         0         0         0         0         0         0
</div><div class="horizontalEllipsis"></div><div class="verticalEllipsis hide"></div></div></div></div><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="9C49FDC1" data-testid="output_9" style="width: 450px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><h2  class = 'S2'><span style=' font-weight: bold;'>Phenotypic phase plane analysis (PhPP)</span></h2><div  class = 'S1'><span>The PhPP is a method for describing in two or three dimensions, how the objective function would change if additional metabolites were given to the model [3]. </span></div><div  class = 'S1'><span>Essentially PhPP performs a </span><span style=' font-family: monospace;'>doubleRobustnessAnalysis()</span><span>, with the difference that shadow prices are retained. The code is as follows-</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >modelphpp = modelalter;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >ATPphppRates = zeros(21);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >i = 0:10</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">for </span><span >j = 0:20</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >        modelphpp = changeRxnBounds(modelphpp, </span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >, -i, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >        modelphpp = changeRxnBounds(modelphpp, </span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >, -j, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >        modelphpp = changeObjective(modelphpp, </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >        FBAphpp = optimizeCbModel(modelphpp, </span><span style="color: rgb(170, 4, 249);">'max'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >        ATPphppRates(i+1,j+1) = FBAphpp.f;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >surfl(ATPphppRates) </span><span style="color: rgb(2, 128, 9);">% 3d plot</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >xlabel(</span><span style="color: rgb(170, 4, 249);">'Flux through EX-glc-D[e]'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >ylabel(</span><span style="color: rgb(170, 4, 249);">'Flux through EX-o2[e]'</span><span >)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre"><span >zlabel(</span><span style="color: rgb(170, 4, 249);">'Objective function'</span><span >)</span></span></div><div  class = 'S4'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="8601E97D" data-testid="output_10" style="width: 450px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><div  class = 'S5'><span>To generate a 2D plot: </span><span style=' font-family: monospace;'>pcolor(ATPphppRates)</span></div><div  class = 'S1'><span>Alternatively, use the function </span><span style=' font-family: monospace;'>phenotypePhasePlane()</span><span>. This function also draws the line of optimality, as well as the shadow prices of the metabolites from the two control reactions. In this case, control reactions are '</span><span style=' font-family: monospace;'>EX_glc_D[e]</span><span>' and '</span><span style=' font-family: monospace;'>EX_o2[e]</span><span>'. The line of optimality signifies the state wherein, the objective function is optimal. In this case it is '</span><span style=' font-family: monospace;'>DM_atp_c_</span><span>'.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >modelphpp = changeObjective (modelphpp, </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >[growthRates, shadowPrices1, shadowPrices2] = phenotypePhasePlane(modelphpp,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >, </span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >);</span></span></div><div  class = 'S4'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="95B686CC" data-testid="output_11" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">generating PhPP</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="652218B2" data-testid="output_12" style="width: 450px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="C02999FA" data-testid="output_13" style="width: 450px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><div  class = 'S5'><span></span></div><h2  class = 'S2'><span>REFERENCES </span></h2><div  class = 'S1'><span>[1] Noronha A., et al. (2017). ReconMap: an interactive visualization of human metabolism. </span><span style=' font-style: italic;'>Bioinformatics</span><span>., 33 (4): 605-607.</span></div><div  class = 'S1'><span>[2] Edwards, J.S. and and Palsson, B. Ø. (2000). Robustness analysis of the Escherichia coli metabolic network. </span><span style=' font-style: italic;'>Biotechnology Progress, </span><span>16(6):927-39.</span></div><div  class = 'S1'><span>[3] Edwards, J.S., Ramakrishna, R. and and Palsson, B. Ø. (2002). Characterizing the metabolic phenotype: A phenotype phase plane analysis. </span><span style=' font-style: italic;'>Biotechnology and Bioengineering</span><span>, 77:27-36.</span></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% Varying Parameters analysis
% *Author(s): Vanja Vlasov, Systems Biochemistry Group, LCSB, University of 
% Luxembourg.*
% 
% *Reviewer(s): Thomas Pfau, Systems Biology Group, LSRU, University of Luxembourg.*
% 
% *Ines Thiele, Molecular Systems Physiology Group, LCSB, University of Luxembourg*
% 
% *Almut Heinken, Molecular Systems Physiology Group, LCSB, University of Luxembourg*
% 
% 
% 
% In this tutorial, we show how computations are performed by varying one or 
% two parameters over a fixed range of numerical values.
%% EQUIPMENT SETUP
% If necessary, initialise the cobra toolbox:

initCobraToolbox(false) % false, as we don't want to update
%% 
% For solving linear programming problems in the analysis, certain solvers are 
% required:

changeCobraSolver ('gurobi', 'all', 1);
%changeCobraSolver ('glpk', 'all', 1);
%% 
% The present tutorial can run with |'glpk'| package, which does not require 
% additional installation and configuration. Although, for the analysis of large 
% models is recommended to use the |'gurobi'| package.
%% PROCEDURE
% Before proceeding with the simulations, the path for the model needs to be 
% set up. In this tutorial, the used model is the generic model of human metabolism, 
% Recon 3 [1]. Therefore, we assume, that the cellular objectives include energy 
% production or optimisation of uptake rates and by-product secretion for various 
% physiological functions of the human body. If Recon 3 is not available, please 
% use Recon 2.

%For Recon3D Change the model
modelFileName = 'Recon2.0model.mat';
modelDirectory = getDistributedModelFolder(modelFileName); %Look up the folder for the distributed Models.
modelFileName= [modelDirectory filesep modelFileName]; % Get the full path. Necessary to be sure, that the right model is loaded
model = readCbModel(modelFileName);

%% 
% If Recon2 is used, the reaction nomenclature needs to be adjusted.

model.rxns(find(ismember(model.rxns,'EX_glc(e)')))={'EX_glc_D[e]'};
model.rxns(find(ismember(model.rxns,'EX_o2(e)')))={'EX_o2[e]'};
%% TROUBLESHOOTING
% If there are multiple energy sources available in the model; Specifying more 
% constraints is necessary. If we do not do that, we will have additional carbon 
% and oxygen energy sources available in the cell and the maximal ATP production. 
% 
% To avoid this issue, all external carbon sources need to be closed.

%Closing the uptake of all energy and oxygen sources
for i=1:length(model.rxns)
    if strncmp(model.rxns{i},'EX_',3)
        model.subSystems{i}='Exchange/demand reaction';
    end
end
idx=strmatch('Exchange/demand reaction', model.subSystems);
c=0;
for i=1:length(idx)
    if model.lb(idx(i))~=0
        c=c+1;
        uptakes{c}=model.rxns{idx(i)};
    end
end

modelalter = model;
modelalter = changeRxnBounds(modelalter, uptakes, 0, 'l');

% The alternative way to do that, in case you were using another large model, 
% that does not contain defined Subsystem is
% to find uptake exchange reactions with following codes:
% [selExc, selUpt] = findExcRxns(model);
% uptakes = model.rxns(selUpt);
% Selecting from the exchange uptake reactions those 
% which contain at least 1 carbon in the metabolites included in the reaction:
% subuptakeModel = extractSubNetwork(model, uptakes);
% hiCarbonRxns = findCarbonRxns(subuptakeModel,1);
% Closing the uptake of all the carbon sources
% modelalter = model;
% modelalter = changeRxnBounds(modelalter, hiCarbonRxns, 0, 'l');
%% Robustness analysis
% Robustness analysis is applied to estimate and visualise how changes in the 
% concentration of an environmental parameter (exchange rate) or internal reaction 
% effect on the objective [2]. If we are interested in varying $$v_{j}$ between 
% two values, i.e., $$v_{j,min}$ and $$v_{j,max}$, we can solve $$l$ optimisation 
% problems:
% 
% $$\begin{array}{ll}\max Z_{k}= c^{T}v \\\text{s.t.} & k =1,...,l,\\& Sv=0,\\\text{fixing} 
% & v_{j}=v_{j,min}+ \frac{(k-1)}{(l-1)}*(v_{j,max}-v_{j,min})\\\text{constraints} 
% & v_{i,min}\leq v_{i} \leq v_{i,max}, i=1,...,n,i \neq j \end{array}$$
% 
% 
% 
% The function |robustnessAnalysis()| is used for this analysis:

% [controlFlux, objFlux] = robustnessAnalysis(model, controlRxn, nPoints,...
%      plotResFlag, objRxn,objType)
%% 
% where inputs are a COBRA model, a reaction that has been analysed and optional 
% inputs:

% INPUTS
% model         COBRA model structure
% controlRxn    Reaction of interest whose value is to be controlled
%
% OPTIONAL INPUTS
% nPoints       Number of points to show on plot (Default = 20)
% plotResFlag   Plot results (Default true)
% objRxn        Objective reaction to be maximized
%               (Default = whatever is defined in model)
% objType       Maximize ('max') or minimize ('min') objective
%               (Default = 'max')
%
% OUTPUTS
% controlFlux   Flux values within the range of the maximum and minimum for
%               a reaction of interest
% objFlux       Optimal values of objective reaction at each control
%               reaction flux value
%% 
% Here, we will investigate how robust the maximal ATP production of the network 
% (i.e., the maximal flux through '|DM_atp_c_'|) is with respect to varying glucose 
% uptake rates and fixed oxygen uptake. 

modelrobust = modelalter;
modelrobust = changeRxnBounds(modelrobust, 'EX_o2[e]', -17, 'b');
AtpRates = zeros(21, 1);
for i = 0:20
    modelrobust = changeRxnBounds(modelrobust, 'EX_glc_D[e]', -i, 'b');
    modelrobust = changeObjective(modelrobust, 'DM_atp_c_');
    FBArobust = optimizeCbModel(modelrobust, 'max');
    AtpRates(i+1) = FBArobust.f;
end
plot (1:21, AtpRates)
xlabel('Flux through EX-glc-D[e]')
ylabel('Objective function')
%% 
% We can also investigate the robustness of the maximal ATP production when 
% the available glucose amount is fixed, while different levels of oxygen are 
% available.

modelrobustoxy = modelalter;
modelrobustoxy = changeRxnBounds(modelrobustoxy, 'EX_glc_D[e]', -20, 'b');
AtpRatesoxy = zeros(21, 1);
for i = 0:20
    modelrobustoxy = changeRxnBounds(modelrobustoxy, 'EX_o2[e]', -i, 'b');
    modelrobustoxy = changeObjective(modelrobustoxy, 'DM_atp_c_');
    FBArobustoxy = optimizeCbModel(modelrobustoxy, 'max');
    AtpRatesoxy(i+1) = FBArobustoxy.f;
end
plot (1:21, AtpRatesoxy)
xlabel('Flux through EX-o2[e]')
ylabel('Objective function')
%% 
% * Double robust analysis*
%% 
% Performs robustness analysis for a pair of reactions of interest and an objective 
% of interest. The double robust analysis is implemented with the function |doubleRobustnessAnalysis()|.

% [controlFlux1, controlFlux2, objFlux] = doubleRobustnessAnalysis(model,...
%  controlRxn1, controlRxn2, nPoints, plotResFlag, objRxn, objType)
%% 
% The inputs are a COBRA model, two reactions for the analysis and optional 
% inputs:

%INPUTS
% model         COBRA model to analyse,
% controlRxn1   The first reaction for the analysis,
% controlRxn2   The second reaction for the analysis;
%
%OPTIONAL INPUTS
% nPoints       The number of flux values per dimension (Default = 20)
% plotResFlag   Indicates whether the result should be plotted (Default = true)
% objRxn        is objective to be used in the analysis (Default = whatever
%               is defined in model)
% objType       Direction of the objective (min or max)
%               (Default = 'max')
%% 
% 

modeldrobustoxy = modelalter;
modeldrobustoxy = changeRxnBounds(modeldrobustoxy, 'EX_glc_D[e]', -20, 'l');
modeldrobustoxy = changeRxnBounds(modeldrobustoxy, 'EX_o2[e]', -17, 'l');
[controlFlux1, controlFlux2, objFlux] = doubleRobustnessAnalysis(modeldrobustoxy,...
    'EX_glc_D[e]', 'EX_o2[e]', 10, 1, 'DM_atp_c_', 'max')
%% *Phenotypic phase plane analysis (PhPP)*
% The PhPP is a method for describing in two or three dimensions, how the objective 
% function would change if additional metabolites were given to the model [3]. 
% 
% Essentially PhPP performs a |doubleRobustnessAnalysis()|, with the difference 
% that shadow prices are retained. The code is as follows-

modelphpp = modelalter;
ATPphppRates = zeros(21);
for i = 0:10
    for j = 0:20
        modelphpp = changeRxnBounds(modelphpp, 'EX_glc_D[e]', -i, 'b');
        modelphpp = changeRxnBounds(modelphpp, 'EX_o2[e]', -j, 'b');
        modelphpp = changeObjective(modelphpp, 'DM_atp_c_');
        FBAphpp = optimizeCbModel(modelphpp, 'max');
        ATPphppRates(i+1,j+1) = FBAphpp.f;
    end
end

surfl(ATPphppRates) % 3d plot
xlabel('Flux through EX-glc-D[e]')
ylabel('Flux through EX-o2[e]')
zlabel('Objective function')
%% 
% To generate a 2D plot: |pcolor(ATPphppRates)|
% 
% Alternatively, use the function |phenotypePhasePlane()|. This function also 
% draws the line of optimality, as well as the shadow prices of the metabolites 
% from the two control reactions. In this case, control reactions are '|EX_glc_D[e]|' 
% and '|EX_o2[e]|'. The line of optimality signifies the state wherein, the objective 
% function is optimal. In this case it is '|DM_atp_c_|'.

modelphpp = changeObjective (modelphpp, 'DM_atp_c_');
[growthRates, shadowPrices1, shadowPrices2] = phenotypePhasePlane(modelphpp,...
    'EX_glc_D[e]', 'EX_o2[e]');
%% 
% 
%% REFERENCES 
% [1] Noronha A., et al. (2017). ReconMap: an interactive visualization of human 
% metabolism. _Bioinformatics_., 33 (4): 605-607.
% 
% [2] Edwards, J.S. and and Palsson, B. Ø. (2000). Robustness analysis of the 
% Escherichia coli metabolic network. _Biotechnology Progress,_ 16(6):927-39.
% 
% [3] Edwards, J.S., Ramakrishna, R. and and Palsson, B. Ø. (2002). Characterizing 
% the metabolic phenotype: A phenotype phase plane analysis. _Biotechnology and 
% Bioengineering_, 77:27-36.
##### SOURCE END #####
-->
</div></body></html>